package ExactMsa;

import java.util.Hashtable;



/*******************************************************************
 * Computing the exact alignment of three sequences with sum-of-pairs-score.
 * Input arguments are three sequences, a score matrix and gap costs.
 *
 * Running time O(n^4*32), where n is length of sequences
 * Space Consumption O(n^4*32)
 *
 * @author Paula C. Tataru, Marcus Schroeder
 *
 *****************************************************************/


public class ExactMsaFourSequences {

	/******************************************************************
	 *  Default constructor
	 *  Computes alignment of hardcoded sequences
	 *****************************************************************/
	public ExactMsaFourSequences(Hashtable<String,Hashtable<String,String>> scoreMatrix,int gapcost) {


		this.scoreMatrix = scoreMatrix;
		this.gapcost = gapcost;

		startAlignment();

	}



	/******************************************************************
	 * Alternative constructor
	 * Computes alignment of three given sequences and a given score matrix
	 *****************************************************************/
	public ExactMsaFourSequences(String s1, String s2, String s3, String s4, Hashtable<String,Hashtable<String,String>> scoreMatrix,int gapcost) {

		this.s1 = s1;
		this.s2 = s2;
		this.s3 = s3;
		this.s4 = s4;


		this.scoreMatrix = scoreMatrix;
		this.gapcost = gapcost;

		startAlignment();

	}




	/******************************************************************
	 * Alternative constructor
	 * Sequences are handed over in an array s
	 *****************************************************************/
	public ExactMsaFourSequences(String[] s, Hashtable<String,Hashtable<String,String>> scoreMatrix,int gapcost) {

		this.s1 = s[0];
		this.s2 = s[1];
		this.s3 = s[2];
		this.s4 = s[3];



		this.scoreMatrix = scoreMatrix;
		this.gapcost = gapcost;

		startAlignment();
	}


	public void startAlignment(){

		double time = System.nanoTime();
		D =  new int[s1.length()+1][s2.length()+1][s3.length()+1][s4.length()+1];

		computeExactAlignment();

		traceback(s1.length(),s2.length(),s3.length(),s4.length());


		double end_time = (System.nanoTime()-time)/1000000000;



		as1 = new StringBuffer(as1).reverse().toString();
		as2 = new StringBuffer(as2).reverse().toString();
		as3 = new StringBuffer(as3).reverse().toString();
		as4 = new StringBuffer(as4).reverse().toString();



		System.out.println("--------------------------------------");
		System.out.println("Results:");
		System.out.println("--------------------------------------");
		System.out.println("Sequence lengths: "+s1.length()+", "+s2.length()+", "+s3.length()+", "+s4.length());
		System.out.println("Taken time in s: "+end_time);
		System.out.println("Score: "+D[s1.length()][s2.length()][s3.length()][s4.length()]);
		System.out.println("--------------------------------------");
		System.out.println();
		System.out.println();
		System.out.println(as1);
		System.out.println(as2);
		System.out.println(as3);
		System.out.println(as4);
		System.out.println();
	}


	/*****************************************************************
	 * GLOBALS
	 *****************************************************************/

	// The four sequences with default values
	private String s1 = "tatggagagaataaaagaactga";
	private String s2 = "atggatgtcaatccgactctact";
	private String s3 = "tccaaaatggaagactttgtgcg";
	private String s4 = "aaaagcaacaaaaatgaaggcaa";

	private String as1 ="";
	private String as2 ="";
	private String as3 ="";
	private String as4 ="";

	/******************************************************************
	 *  Gapcosts
	 *****************************************************************/
	private int gapcost = -5;


	/******************************************************************
	 *  The Score Matrix
	 *****************************************************************/
	private Hashtable<String,Hashtable<String,String>> scoreMatrix;







	/******************************************************************
	 *  Alignment Matrix D
	 *
	 *  Has four dimensions with size [n,m,o,p], where
	 *  			n is length of sequence 1,
	 *  			m is length of sequence 2,
	 *  			o is length of sequence 3,
	 *  			p is legnth of sequence 4
	 *
	 *  Is intialized in complete size (n*m*o*p) on class initialization.
	 *****************************************************************/
	private int[][][][] D;



	/******************************************************************
	 * Returns the score for two compared characters a and b
	 * Values stems from a score matrix
	 *****************************************************************/
	private int matchScore(char a,char b){

		String score = ((Hashtable<String,String>) scoreMatrix.get(""+a)).get(""+b).toString();

		return Integer.parseInt(score);

		//return Float.parseFloat(score);

	}


	/******************************************************************
	 * Compute exact alignment of three sequences.
	 *****************************************************************/
	private void computeExactAlignment(){

		int n = s1.length();
		int o = s2.length();
		int p = s3.length();
		int q = s4.length();



		/******************************************************
		 *  INTIALIZING THE ALIGNMENT MATRIX
		 ******************************************************/
		D[0][0][0][0] = 0;


		// D_s1_s2_s3

		ExactMsaThreeSequences d_sub = new ExactMsaThreeSequences(s1,s2,s3,scoreMatrix,gapcost);

		for(int i = 0; i <= n; i++){
			for(int j = 0; j <= o; j++){
				for(int k = 0; k <= p; k++){

					D[i][j][k][0] = d_sub.getValueInAlignmentMatrix(i, j, k) + (i+j+k)*gapcost;

				}
			}
		}

		// D_s1_s3_s4

		d_sub = new ExactMsaThreeSequences(s1,s3,s4,scoreMatrix,gapcost);

		for(int i = 0; i <= n; i++){
			for(int k = 0; k <= p; k++){
				for(int m = 0; m <= q; m++){

					D[i][0][k][m] = d_sub.getValueInAlignmentMatrix(i, k, m) + (i+k+m)*gapcost;

				}
			}
		}

		// D_s1_s2_s4
		d_sub = new ExactMsaThreeSequences(s1,s2,s4,scoreMatrix,gapcost);

		for(int i = 0; i <= n; i++){
			for(int j = 0; j <= o; j++){
				for(int m = 0; m <= q; m++){

					D[i][j][0][m] = d_sub.getValueInAlignmentMatrix(i, j, m) + (i+j+m)*gapcost;

				}
			}
		}

		// D_s2_s3_s4
		d_sub = new ExactMsaThreeSequences(s2,s3,s4,scoreMatrix,gapcost);

		for(int j = 0; j <= o; j++){
			for(int k = 0; k <= p; k++){
				for(int m = 0; m <= q; m++){

					D[0][j][k][m] = d_sub.getValueInAlignmentMatrix(j,k, m) + (j+k+m)*gapcost;

				}
			}
		}


		/******************************************************
		 *  FILLING THE TABLE
		 ******************************************************/

		for(int i = 1; i <= n; i++){
			for(int j = 1; j <=  o; j++){
				for(int k = 1; k <= p; k++){
					for(int m = 1; m <= q; m++){


						int c_ij = matchScore(s1.charAt(i-1), s2.charAt(j-1));
						int c_ik = matchScore(s1.charAt(i-1), s3.charAt(k-1));
						int c_im = matchScore(s1.charAt(i-1), s4.charAt(m-1));
						int c_jk = matchScore(s2.charAt(j-1), s3.charAt(k-1));
						int c_jm = matchScore(s2.charAt(j-1), s4.charAt(m-1));
						int c_km = matchScore(s3.charAt(k-1), s4.charAt(m-1));

						int d1 = D[i-1][j-1][k-1][m-1] + c_ij + c_ik + c_im + c_jk + c_jm + c_km;

						int d2 = D[i-1][j-1][k-1][m] + c_ij + c_ik + c_jk + 3*gapcost;
						int d3 = D[i-1][j-1][k][m-1] + c_ij + c_im + c_jm + 3*gapcost;
						int d4 = D[i-1][j][k-1][m-1] + c_ik + c_im + c_km + 3*gapcost;
						int d5 = D[i][j-1][k-1][m-1] + c_jk + c_jm + c_km + 3*gapcost;


						int d6  = D[i-1][j-1][k][m] + c_ij + 4*gapcost;
						int d7  = D[i-1][j][k][m-1] + c_im + 4*gapcost;
						int d8  = D[i-1][j][k-1][m] + c_ik + 4*gapcost;
						int d9  = D[i][j-1][k][m-1] + c_jm + 4*gapcost;
						int d10 = D[i][j-1][k-1][m] + c_jk + 4*gapcost;
						int d11 = D[i][j][k-1][m-1] + c_km + 4*gapcost;

						int d12 = D[i][j][k][m-1] + 3*gapcost;
						int d13 = D[i][j][k-1][m] + 3*gapcost;
						int d14 = D[i][j-1][k][m] + 3*gapcost;
						int d15 = D[i-1][j][k][m] + 3*gapcost;


						// 	FINDING THE MINIMUM ENTRY
						D[i][j][k][m] = Math.min(
											Math.min(
													Math.min(
															Math.min(d1, d2),
															Math.min(d3, d4)
													),
													Math.min(
															Math.min(d5,d6),
															Math.min(d7,d8)
													)
											),
											Math.min(
													Math.min(
															Math.min(d9, d10),
															Math.min(d11, d12)
													),
													Math.min(
															Math.min(d13,d14),
															d15
													)
											)
									);
					}
				}

			}

		}
	}


	private void traceback(int i, int j, int k, int m){


            if(i == 0 && j == 0 && k == 0 && m == 0){

		}
		else{
			int c_ij = 0;
			int c_ik = 0;
			int c_im = 0;
			int c_jk = 0;
			int c_jm = 0;
			int c_km = 0;

			if(i>0 && j >0) c_ij = matchScore(s1.charAt(i-1), s2.charAt(j-1));
			if(i>0 && k >0) c_ik = matchScore(s1.charAt(i-1), s3.charAt(k-1));
			if(i>0 && m >0) c_im = matchScore(s1.charAt(i-1), s4.charAt(m-1));
			if(j>0 && k >0) c_jk = matchScore(s2.charAt(j-1), s3.charAt(k-1));
			if(j>0 && m >0) c_jm = matchScore(s2.charAt(j-1), s4.charAt(m-1));
			if(k>0 && m >0) c_km = matchScore(s3.charAt(k-1), s4.charAt(m-1));

			if((i>0 && j>0 && k>0 && m>0) && D[i][j][k][m] == D[i-1][j-1][k-1][m-1] + c_ij + c_ik + c_im + c_jk + c_jm + c_km){

					as1 = as1+s1.charAt(i-1);
					as2 = as2+s2.charAt(j-1);
					as3 = as3+s3.charAt(k-1);
					as4 = as4+s4.charAt(m-1);

					traceback(i-1,j-1,k-1,m-1);

			}
			else if((i>0 && j>0 && k>0) && D[i][j][k][m] == D[i-1][j-1][k-1][m] + c_ij + c_ik + c_jk + 3*gapcost){

					as1 = as1+s1.charAt(i-1);
					as2 = as2+s2.charAt(j-1);
					as3 = as3+s3.charAt(k-1);
					as4 = as4+"-";

					traceback(i-1,j-1,k-1,m);
			}
			else if((i>0 && j>0 && m>0) && D[i][j][k][m] == D[i-1][j-1][k][m-1] + c_ij + c_im + c_jm + 3*gapcost){

				as1 = as1+s1.charAt(i-1);
				as2 = as2+s2.charAt(j-1);
				as3 = as3+"-";
				as4 = as4+s4.charAt(m-1);

				traceback(i-1,j-1,k,m-1);
			}
			else if((i>0 && k>0 && m>0) && D[i][j][k][m] == D[i-1][j][k-1][m-1] + c_ik + c_im + c_km + 3*gapcost){

				as1 = as1+s1.charAt(i-1);
				as2 = as2+"-";
				as3 = as3+s3.charAt(k-1);
				as4 = as4+s4.charAt(m-1);

				traceback(i-1,j,k-1,m-1);
			}
			else if((j>0 && k>0 && m>0) && D[i][j][k][m] == D[i][j-1][k-1][m-1] + c_jk + c_jm + c_km + 3*gapcost){

				as1 = as1+"-";
				as2 = as2+s2.charAt(j-1);
				as3 = as3+s3.charAt(k-1);
				as4 = as4+s4.charAt(m-1);

				traceback(i,j-1,k-1,m-1);
			}
			else if((i>0 && j>0) && D[i][j][k][m] == D[i-1][j-1][k][m] + c_ij + 4*gapcost){

				as1 = as1+s1.charAt(i-1);
				as2 = as2+s2.charAt(j-1);
				as3 = as3+"-";
				as4 = as4+"-";

				traceback(i-1,j-1,k,m);
			}
			else if((i>0 && k>0) && D[i][j][k][m] == D[i-1][j][k-1][m] + c_ik + 4*gapcost){

				as1 = as1+s1.charAt(i-1);
				as2 = as2+"-";
				as3 = as3+s3.charAt(k-1);
				as4 = as4+"-";

				traceback(i-1,j,k-1,m);
			}
			else if((i>0 && m>0) && D[i][j][k][m] == D[i-1][j][k][m-1] + c_im + 4*gapcost){

				as1 = as1+s1.charAt(i-1);
				as2 = as2+"-";
				as3 = as3+"-";
				as4 = as4+s4.charAt(m-1);

				traceback(i-1,j,k,m-1);
			}
			else if((j>0 && k>0) && D[i][j][k][m] == D[i][j-1][k-1][m] + c_jk + 4*gapcost){

				as1 = as1+"-";
				as2 = as2+s2.charAt(j-1);
				as3 = as3+s3.charAt(k-1);
				as4 = as4+"-";

				traceback(i,j-1,k-1,m);
			}
			else if((j>0 && m>0) && D[i][j][k][m] == D[i][j-1][k][m-1] + c_jm + 4*gapcost){

				as1 = as1+"-";
				as2 = as2+s2.charAt(j-1);
				as3 = as3+"-";
				as4 = as4+s4.charAt(m-1);

				traceback(i,j-1,k,m-1);
			}
			else if((k>0 && m>0) && D[i][j][k][m] == D[i][j][k-1][m-1] + c_km + 4*gapcost){

				as1 = as1+"-";
				as2 = as2+"-";
				as3 = as3+s3.charAt(k-1);
				as4 = as4+s4.charAt(m-1);

				traceback(i,j,k-1,m-1);
			}
			else if(i>0 && D[i][j][k][m] == D[i-1][j][k][m] + 3*gapcost){

				as1 = as1+s1.charAt(i-1);
				as2 = as2+"-";
				as3 = as3+"-";
				as4 = as4+"-";

				traceback(i-1,j,k,m);
			}
			else if(j>0 && D[i][j][k][m] == D[i][j-1][k][m] + 3*gapcost){

				as1 = as1+"-";
				as2 = as2+s2.charAt(j-1);
				as3 = as3+"-";
				as4 = as4+"-";

				traceback(i,j-1,k,m);
			}
			else if(k>0 && D[i][j][k][m] == D[i][j][k-1][m] + 3*gapcost){

				as1 = as1+"-";
				as2 = as2+"-";
				as3 = as3+s3.charAt(k-1);
				as4 = as4+"-";

				traceback(i,j,k-1,m);
			}
			else if(m>0 && D[i][j][k][m] == D[i][j][k][m-1] + 3*gapcost){

				as1 = as1+"-";
				as2 = as2+"-";
				as3 = as3+"-";
				as4 = as4+s4.charAt(m-1);

				traceback(i,j,k,m-1);
			}




		}


	}

}
